function  dreamMaster(dreamPar)
  global totalNumberOfNodes; global DIE_TAG;
   
    %initDreamPar;
    
    nElemOutput = ceil(dreamPar.nModelEvalsMax/dreamPar.nSeq);
    
    iteration = 0;

    outliers=[];
    logpHistory = nan(floor(dreamPar.nModelEvalsMax/dreamPar.nSeq),dreamPar.nSeq);
    deltaTot = zeros(1,dreamPar.nCrossoverValues);
    pCrossoverHistory = nan(floor(dreamPar.nModelEvalsMax/dreamPar.nSeq),dreamPar.nCrossoverValues +1);
    critGelRub = nan(nElemOutput, dreamPar.nOptPars+1);
    acceptanceRate = nan(nElemOutput,2);
    acceptanceRate(1,1:2) = [dreamPar.nSeq - 1];

    if dreamPar.crossoverValuesTuning
        % Calculate multinomial probabilities of each of the nCrossoverValues 
        pCrossover = (1/dreamPar.nCrossoverValues) * ones(1,dreamPar.nCrossoverValues);
        % Calculate the actual allCrossoverValues values based on p
        [allCrossoverValues] = generateCrossoverValues(dreamPar,pCrossover); 
        crossoverCount = zeros(1,dreamPar.nCrossoverValues);
    else
        pCrossover = 1/dreamPar.nCrossoverValues; 
        crossoverCount = [];
        [allCrossoverValues] = pCrossover * ones(dreamPar.nSeq,dreamPar.nOffspringPerSeq); 
        crossoverCount = dreamPar.nSeq * dreamPar.nOffspringPerSeq;
    end;
    
    for i = 1:dreamPar.nDiffEvolPairs,
        jumpRateTable(:,i) = 2.38./sqrt(2 * i * [1:dreamPar.nOptPars]'); 
    end;
    
    pCrossoverHistory(1,1:size(pCrossover,2)+1) = [1 pCrossover]; 


    switch dreamPar.samplingMethod  
        case 'covariance'
            randomDraw = covRandomDraw(dreamPar);
            evalResults = [repmat(NaN,[dreamPar.nSamples,1]),randomDraw,...
                           repmat([NaN,NaN,false],[dreamPar.nSamples,1])];
            nModelEvals = 0;               
            [evalResults,nModelEvals] = calcObjectiveScoreParallel(dreamPar,evalResults,totalNumberOfNodes, nModelEvals,iteration);

            % partition 'evalResults' into sequences and initialize the sequences:
            sequences = toSequences(dreamPar,evalResults);
        case 'uniform'
            % Uniform sampling of points in parameter space:
            randomDraw = stratranddraw(dreamPar);
            evalResults = [repmat(NaN,[dreamPar.nSamples,1]),randomDraw,...
                           repmat([NaN,NaN,false],[dreamPar.nSamples,1])];
            nModelEvals = 0;               
            [evalResults,nModelEvals] = calcObjectiveScoreParallel(dreamPar,evalResults,totalNumberOfNodes, nModelEvals,iteration);

            % partition 'evalResults' into sequences and initialize the sequences:
            sequences = toSequences(dreamPar,evalResults);

        case 'continue'
            % continue with an earlier run
            reply = input([char(10),'Please provide the filename of the *.mat ',...
                           'file containing',char(10),'variables ',...
                            char(39),'critGelRub',char(39),', ',...
                            char(39),'evalResults',char(39),', and ',...
                            char(39),'sequences',char(39),'.',char(10),'>> ']);
            if exist(reply,'file')
                eval(['load(',char(39),reply,char(39),',',...
                    char(39),'critGelRub',char(39),',',...
                    char(39),'evalResults',char(39),',',...
                    char(39),'sequences',char(39),')'])
            else
                error('No such file exists on the path.')
            end

            nRecords = max(evalResults(:,dreamPar.iterCol));

            % partition 'evalResults' into sequences:
            newSequences = toSequences(dreamPar,evalResults);
            dreamPar.converged = abortoptim(critGelRub,dreamPar)    
            % initialize the sequences:
            if isempty(newSequences)
                sequences = newSequences;
                nModelEvals = max(evalResults(:,dreamPar.iterCol));
            else
                nTrownAway = mod(nRecords-dreamPar.nSamples,dreamPar.nOffspring);
                nModelEvals = nRecords-nTrownAway;
                if dreamPar.verboseOutput            
                    disp(['Throwing away last ',num2str(nTrownAway), ' model evaluations.'])
                    disp(['dream run continues at ',num2str(nModelEvals+1), ' model evaluations.'])    
                end
            end

    end
    
     % Save history log density of individual chains
   % nRecords = max(evalResults(:,dreamPar.iterCol));
   
  
    logpHistory(1,1:nModelEvals+1) = [1 evalResults(:,dreamPar.logPCol)'];
    critGelRub(1,:) = [1,gelmanrubin(dreamPar,sequences)];
%     evalResults
    lastPointsFromEverySeq = getLastPointsFromEverySeq(dreamPar,sequences);
    
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %    
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %    
    % % % % % %                                             % % % % % %    
    % % % % % %         DREAM  START OF MAIN LOOP         % % % % % %    
    % % % % % %                                             % % % % % %    
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %    
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %   
     
    [sequences,evalResults] = preAllocate(dreamPar,sequences,evalResults);
%     return;
    while continueDream(dreamPar,nModelEvals)
  
            for iOffspring = 1:dreamPar.nOffspringPerSeq
                
                iteration = iteration + 1;

                % propose offspring 
                propChild = generateOffspring(dreamPar,lastPointsFromEverySeq, jumpRateTable, allCrossoverValues(:,iOffspring));

                % calculate objective for unevaluated rows of evalResults:
                [propChild,nModelEvals] = calcObjectiveScoreParallel(dreamPar,propChild, totalNumberOfNodes,nModelEvals, iteration);

                % evolve the sequences
                [acceptedChild,alpha,acceptedSequences] = evolve(dreamPar,lastPointsFromEverySeq,propChild);
               
                rejectedSeqIdx = find(acceptedSequences==0);

                % Only apply delayed rejection to the chains that were
                % rejected in the previous phase 
                if dreamPar.delayedRejection && numel(rejectedSeqIdx) > 0
                    
                    propChildDelayed = generateOffspringDelayedStep(dreamPar,lastPointsFromEverySeq);
                    
                    %nModelEvals = nModelEvals - dreamPar.nSeq; %size(propChildDelayed,1);
                    
                    [propChildDelayed,nModelEvals] = calcObjectiveScoreParallel(dreamPar,propChildDelayed(rejectedSeqIdx,:),totalNumberOfNodes,nModelEvals, iteration);
                    
                    [acceptedChildDelayed,alphaDelayed,acceptedSequencesDelayed] = evolve(dreamPar,propChildDelayed, propChild(rejectedSeqIdx,:));
                    
                    [acceptedChild,acceptedSequencesDelayed] = evolveDelayedStep(dreamPar,lastPointsFromEverySeq, propChildDelayed, propChild, acceptedChild,rejectedSeqIdx,alpha,alphaDelayed);
                    
                end;

                % add the accepted child to the 'evalResults' array:
                 evalResults(iteration*dreamPar.nSeq+1:(iteration+1)*dreamPar.nSeq, :) = acceptedChild;

                
                % add the accepted child to the 'sequences' array:
               
                [sequences,wasLastPointStored] = updateSequences(dreamPar, sequences, acceptedChild, iOffspring, iteration+1);
                if(dreamPar.reducedSampleCollection) 
                        updateIndex = iteration/dreamPar.reducedSampleFrequency;
                    else updateIndex = iteration;
                end;
                
                if dreamPar.crossoverValuesTuning
                    [deltaTot] = updateDelta(dreamPar,deltaTot,acceptedChild, lastPointsFromEverySeq,allCrossoverValues(:,iOffspring));
                end;
                
                lastPointsFromEverySeq = acceptedChild;
               
                logpHistory(iteration,:) = [iteration * dreamPar.nSeq acceptedChild(:,dreamPar.logPCol)']; %iteration*dreamPar.nSeq 
                if(dreamPar.delayedRejection) && numel(rejectedSeqIdx) > 0
                    acceptanceRate(iteration,1:2) = [iteration*dreamPar.nSeq 100 * (sum(acceptedSequences) + sum(acceptedSequencesDelayed))/(dreamPar.nSeq + size(rejectedSeqIdx,1))];
                else
                    acceptanceRate(iteration,1:2) = [iteration*dreamPar.nSeq 100 * sum(acceptedSequences)/dreamPar.nSeq ];
                end;
                
            end

           


        % Check whether to update individual pCrossover values
        if (iteration*dreamPar.nSeq <= 0.1 * dreamPar.nModelEvalsMax)
            if dreamPar.crossoverValuesTuning
                [pCrossover,crossoverCount] = updatePCrossover(dreamPar,allCrossoverValues,deltaTot,crossoverCount); 
            end;
        else
            % See whether there are any outlier chains, and replace them to
            % current best value of lastPointsFromEverySeq
            [lastPointsFromEverySeq,logpHistory,outliers] = removeOutlierChains(lastPointsFromEverySeq,...
                                                                     logpHistory,iteration*dreamPar.nSeq,outliers,dreamPar);
            if wasLastPointStored
                for seq = 1:dreamPar.nSeq 
                    sequences(updateIndex,:,seq) = lastPointsFromEverySeq(seq, :);
                end
            end;
        end;
        % Store Probability of individual crossover values
        pCrossoverHistory(floor(iteration / dreamPar.nOffspringPerSeq)+1,:) = [iteration pCrossover];

        if dreamPar.crossoverValuesTuning
            [allCrossoverValues] = generateCrossoverValues(dreamPar,pCrossover); 
        else
            [allCrossoverValues] = pCrossover * ones(dreamPar.nSeq,dreamPar.nOffspringPerSeq); 
        end;

       %determine the Gelman-Rubin scale reduction (convergence) statistic:
        critGelRub(floor(iteration / dreamPar.nOffspringPerSeq)+1,:) = [iteration,gelmanrubin(dreamPar,sequences(1:updateIndex,:,:))];
            
        % determine whether all parameters meet the Gelman-Rubin scale
        % reduction criterion
        dreamPar.converged = abortoptim(critGelRub,dreamPar);
  

    end
    
    
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %    
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %    
    % % % % % %                                             % % % % % %    
    % % % % % %          DREAM END OF MAIN LOOP           % % % % % %    
    % % % % % %                                             % % % % % %    
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %    
    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %    

    
   firstNan = find(isnan(evalResults(:,1)),1,'first');
   evalResults = evalResults(1:firstNan-1,:);
   firstNan = find(isnan(acceptanceRate(:,1)),1,'first');
   acceptanceRate = acceptanceRate(1:firstNan-1,:);
   firstNan = find(isnan(sequences(:,1,1)),1,'first');
   sequences = sequences(1:firstNan-1,:,:);
   firstNan = find(isnan(critGelRub(:,1)),1,'first');
   critGelRub = critGelRub(1:firstNan-1,:);
    
    for worker = 1:totalNumberOfNodes-1
        package = [];
        MPI_Send(package,worker, DIE_TAG, MPI_COMM_WORLD);
    end


    save('outputDREAM.mat', 'evalResults','critGelRub','outliers','acceptanceRate', 'pCrossoverHistory', 'sequences','dreamPar');



